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We consider abstract operator equations Fu — y, where F is a compact 
linear operator between Hilbert spaces U and V, which are function spaces 
on closed, finite dimensional Riemannian manifolds, respectively. This 
setting is of interest in numerous applications such as Computer Vision 
and non-destructive evaluation. 

In this work, we study the approximation of the solution of the ill-posed 
operator equation with Tikhonov type regularization methods. We prove 
well-posedness, stability, convergence, and convergence rates of the regu- 
larization methods. Moreover, we study in detail the numerical analysis 
and the numerical implementation. Finally, we provide for three different 
inverse problems numerical experiments. 

Key words: Inverse problems, variational regularization on Rieman- 
nian manifolds, functions of bounded variation 

1 Introduction 

The problem of solving linear inverse and ill-posed problems has a long 
tradition in engineering (see [16]). Several strategies have been proposed in the 
literature to solve such problems approximatively in a stable manner. 

However, in most applications the data are assumed to be functions, which 
are defined on a subset of an Euclidean space. In this paper the focus is on 
imaging problems, where the data are functions on closed, finite dimensional 
Riemannian manifolds. Such problems appear in Computer Vision and non- 
destructive evaluation, to name but a few (cf. Section [3j). 

In this paper we take an abstract point of view and formulate the ill-posed 
imaging problem as the solution of an operator equation 
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Abstract 



Fu = y 



(1) 



Here F describes the physics of image formation, and y denotes the ideal mea- 
surement data, which contains neither noise 5 nor modeling errors. The operator 
F : U — > V is a compact linear operator between Hilbert spaces of functions 
defined on closed, finite dimensional Riemannian manifolds, respectively. Con- 
sequently, the inverse operator is unbounded and the solution of ([l]) is ill-posed. 

In practice, ideal data are not available, but rather some approximation y 5 . 
These perturbations, in general, do not allow for a direct stable inversion of 
F. To provide a stable numerical solution, Tikhonov type regularization is an 
adequate choice (see [T5l [22l [TTJ [27l E6]). This method consisting in calculating 
a minimizer u s a of the functional 

T atyS (u):=^\\Fu-y s \\ 2 + an(u), (2) 

which approximates the solution of ([I]). Here, typically, TZ : U — > [0, +oo] is a 
proper, convex regularization functional. The parameter a controls the trade 
off between the quality of approximation of Fu and y s as well as the stability 
of the minimizer. The choice of the regularizing functional TZ is essential and 
is selected according to problem specifications. Typical choices of TZ, which 
are also considered in the paper, are Sobolev space (semi-)norms and the total 
variation semi-norm on manifolds. In this paper, we do not consider more 
general settings of non-convex regularization functionals, as it has been done in 
the Euclidean setting for instance in [351 03] • 

In this paper we consider three different applications of regularization meth- 
ods on manifolds, which are denoising, deblurring, and an inverse problem from 
non-destructive evaluation, which has been studied recently in [20 . 

In the following we summarize some related work: Diffusion filtering on sur- 
faces has been used successfully for denoising [5J [5] , which can be considered 



a particular inverse problem (see Section 3.2 1. Even more multi-scale decom- 
position of data on manifolds can be used for denoising [U [12] . The numerical 
analysis and implementation of this paper is related to the work for discretiza- 
tion of partial differential equations on manifolds, in particular discretization 
of the Laplace-Beltrami operator on the manifold M, — A_a^. Here, in partic- 
ular, we refer to pioneering work of Dzuik 9 on surface finite elements. The 
estimates there were generalized [B] by considering adaptive finite elements. 
Subsequently, several parabolic diffusion equations, (THl E] (isotropic) and [5] 
(anisotropic), were developed for manifold valued data. This topic should not 
be confused with the topic of the paper, where the domain of the functions is a 
manifold, where the context of the other papers is that the functions range in a 
manifold (see [21]). 

The paper is organized as follows: In Section [2] we prove well-posedness of 
variational regularization on closed, finite dimensional Riemannian manifolds. 
Also, convergence rates with respect to Bregman distances are obtained in the 
convex regularization setting, under a standard source condition. Section [3] is 
concerned with numerical minimization of the discrete Tikhonov functional - 
this is most probably the most important contribution of this paper. We pro- 
vide a consistent discretization of convex Tikhonov functionals and formulate 
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them in a purely matrix analysis fashion. As a byproduct this approach pro- 
vides a consistent discretization of some nonlinear partial differential operators. 
Moreover, the consistent discretization is the basis to solve inverse problems 
in a stable way. Section |3.2| provides numerical experiments for three different 
applications. Finally, in Section [5j we provide the basic notions on differential 
geometry and non-linear analysis on manifolds and provide some embedding 
results for Sobolev space and the space of functions of finite total variation. 

2 Analysis of variational regularization for func- 
tions on Riemannian manifolds 

In this section, we state an analysis of variational regularization methods 
for solving the ill-posed operator Equation ([I]) for functions on manifolds. Well 
definedness, stability, convergence, and convergence rate are proven along the 
lines of [55] - the manifold setting does not further complicate the analysis, 
and thus is omitted. However, the results are formulated below for the sake of 
completeness and fixation of the notation: 

Assumption 2.1. (AX) U and V are Hilbert spaces andru ,t\? denote the weak 
topologies, respectively. 

(A2) The functional TZ : U — > [0, +oo] is convex and sequentially lower semi- 
continuous with respect to Tjj ■ 

(A3) T> := T>(F) n T>(1Z) ^ (which in particular implies that TZ is proper). 

(AA) For every a > and C > 0, the lower level set of the Tikhonov functional 

\eve\ c (T a ^) := {u e U : T^ yS (u) < C} 

are sequentially pre-compact with respect to tjj- 

(A5) For every a > and C > 0, the set \evelc(T a y s) is sequentially closed 
with respect to tjj and the restriction of F to levelc(T a y s) is sequentially 
continuous with respect to tjj and Ty . 

The results from |26j imply then: 

Theorem 2.2. Let Assumption \2. 1\ hold. Then, 

• there exists a minimizer of T a y s for every a > and y s 6 V. 

• Let a > 0. Then, for every sequence yt — > y s let us denote 

u k G argminT Qi y fc , k e N, 

then (ufc) has a convergent subsequence. Every convergent subsequence of 
(uk) converges to a minimizer ofT a y s. 
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• Assume that Equation has a solution in T> . Moreover, assume that a 
function a : (0, oo) — > (0,oo) satisfies 

5 2 

a(5) -> and -— - -> for 5 ->• 0. 
a(d) 

Lei i/ie sequence 5 k of positive numbers converging to and assume that 
the data y k := y Sk , a k := a(5 k ) satisfies \\y - y k \ < S k . 

Then, u k G argminT^^)^ has a convergent subsequence and every limit 
is a solution of Equation Q). 

For obtaining qualitative estimates for the convergence of a Tikhonov reg- 
ularized solution to a minimum norm solution, some additional assumptions, 
such as the so-called source condition, are needed. 



Proposition 2.3 (Convergence rates). Let Assumption 2.1 hold. Assume there 
exists a 7Z minimizing solution G D(F) n V(1Z) of Equation and an 
element 

f G Ran(F*) fl^u 1 ). (3) 
Then, with the parameter choice a ~ 5, we have 

Ds(u 5 a ,vl) = 0(5) and \\Fu 5 a -y s \\= 0(5) , 
where denotes the Bregman distance, which is defined as follows: 

D^(u, u) := TZ(u) — TZ(u) — (£, u — u) , u,u € U . 
Here (£, u — u) denotes the inner product on the Hilbert space U . 

Proposition |2.3| applies for instance to total variation minimization 

T a>v i (u) := \ I (Fu -y 5 ) 2 + a \D M u\ (M) , 

1 J M 



where 

\D M u\ (M) 



= sup |y w!yv m X dv(g) : X G C°°(M,R n ) and |JC| L -, (A<iH «) < lj 



(4) 



denoted the total variation of u on the manifold and V» denotes the covariant 
derivative. We choose the space U, V = L 2 (A4). Moreover, we assume that 
F is continuous on L 2 (M) with V(F) n V(K) ^ 0. K(u) = \D M u\ (M) is 
the total variation semi-norm. The verification of Assumption |2.l| is similar to 
the Euclidean setting, and thus omitted. However, the verification requires the 



Meyer-Serrin Theorem |5.2| and the Compactness Theorem 5.3 for functions of 
Bounded Variation BV(Ai) on manifolds. Using both theorems allows to shows 
that the Poincare inequality holds and from this follows that T ay s is coercive 
(A. 4). The Compactness Theorem is applied to verify (A. 5). 

Interpretation of the source condition and the convergence rates have been 
given in [26] for function defined on subsets of M", but are valid in the manifold 
setting in a completely analogous manner. 
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3 Numerical results 



In this section we discuss the implementation of variational regularization 
method for functions defined on manifolds. Afterwards three inverse problems 
and numerical experiments are considered. The three applications are denoising, 
deblurring, and an inverse problem for the Funk-Radon transform. 

Now we discuss the numerical minimization of the discrctized Tikhonov func- 
tional. 

We assume that the closed Riemannian manifold M is approximated by a 
polyhedron M represented as M = (V, T), with vertices V = {vi, . . . , vk} G K 3 
and triangles T = {Xi, . . . , Tj,} cVxVxV. The three components of a vertex 
Vk are denoted by vi with j = 1, . . . , 3. Each triangle Tj € T is defined by the 
set of indices tj = {ii, *2, *3}, of the vertices {vi 1 , Vi 2 , Uj 3 }, which are assumed to 
be counter-clockwise oriented. In this section we only deal with the manifold M , 
and assume that M is a sufficiently good approximation to M. which justifies 
an identification. Consequently, also the metric g and the surface measure on 
M, v{g), are also identified. 

The polyhedral surfaces used in the numerical experiments below have been 
taken from the database [23] • Each surface is closed, of genus zero, and consists 
of approximately 25000 vertices. For genus zero surfaces a natural parametriza- 
tion is the sphere. Following [13] imaging testdata y s on the manifold M is 
generated by mapping a given function with planar domain Q C K 2 onto M by 
making use of the spherical parametrization. 

Polyhedral representation 

Each triangle Tj is parameterized with respect to its vertices Vk, Vj, and vi 
by using barycentric coordinates 

*k,i(7) = v k + (i(vj - Vk) + C, 2 (vi - Vk) , 

where 

T = {7 = (Ci, Ca) : Ci € [0, 1] and ( 2 € [0, 1 - Ci]} • 

We approximate the minimizer of the Tikhonov functional from ^ by the 
minimizer of T a y s on the finite dimensional space of piecewise linear functions 
on the polyhedron M: For k £ {1, . . . , K} let (<Pk) be the function, which is 
continuous on M , linear on each triangle T%,i = \,...,L, and satisfies <fk{vk) — 
1 and tpk{v s ) = if s k. On each triangle Tj we have exactly three such 
functions (for every vertex). The set of piecewise linear functions is the linear 
span of the functions (fk- 

PL(M) := • 

From the definition of ipk and Xk,i it follows that 

x = Vk^Pki^) for every x G M. (5) 
fc 
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and 

iPfc(xfc,i(7)) = 1 - Ci - C2 for every 7 <E T . (6) 

Minimization of the Tikhonov functional is performed for u € PL(M) and we 
assume that the data y s e PL(M) too. Thus the functions u over which we 
minimize and the data y s can be expressed via there series expansion: 



M W = ^2u k f k (x) and y s (x) = ^y^ fe (x) . 



(7) 



The vectors of coefficients are denoted in boldface by u := («fe), y s ■= (yl), 
respectively. The Jacobian of the parametrization of the manifold M is the 
matrix 

/ Ji ... \ 



Jo 



V 



p 2Lx3L 



(8) 



with L blocks 



Ji 







9<Pi 


Ci 




Ci 






dipi 


. C2 


C2 


C2 . 





' -1 1 


" 




-1 


1 



Each submatrix J^ is the Jacobian of the parameterizations in the triangle Tj. 
All vertices in M are put in a block diagonal matrix V <E M 3Lx3i with L blocks 



Each submatrix Vj stores the vertices from triangle Tj ordered accordingly to 
the basis functions Lp k . The metric tensor G € R 2Lx2L on M is block diagonal 
matrix with the L diagonal blocks 



(V.jf ) T (V,Jf) € 



d2x2 



Again, G, is the metric tensor in a given triangle Tj. Let Ai denote the area of 
the triangle Tj, then the volume of the metric tensor satisfies: 

- V\dct(G t )\ = 2A, , 

Therefore the surface measure dv(g) can be expressed in barycentric coordinates 
and the relation reads as follows: 

du(g) = 2A t dl- 

Let V G M. 3LxK a matrix which encodes the connectivity of the manifold M. 
That is 



V 



ki 



1 if i e t fe / where fc = 3(fc' - 1) + j, j = 1, 2, 3 
else. 



G 



Vi £ M. 3xK is a linear mapping assigning each triangle T, the indices of the 
three vertices. Accordingly, the covariant derivative of u e PL(M), Vm«, on 
the triangle Ti is a constant vector and is given by 



V,3f G^ViU G 



c3xl 



(9) 



The matrix 



/ Zi 





V o 





z 2 











p 3LxL 



consists of the gradient vectors of u on each triangle of M. The matrix Z T Z is a 
positive semi-definite diagonal matrix. Thus (Z T Z) p / 2 is the matrix consisting 
of the p/2 powers of diagonal entries. 



3.1 Discretization of the Tikhonov Functional 

In the following we consider minimization of the discrete Tikhonov functional 
with functions defined on PL(M). The goal is to express the fit-to-data term 
and the rcgularization functional in dependence of the vector u. 

In all our test cases we have that F : U — > V, where U and V are function 
spaces defined on the same closed, finite dimensional Riemannian manifold M. 
We assume that Fu can be approximated by a piecewise linear function on the 
polyhedron M (note that here both Fu and M. are approximated). 

The linear operator F may not necessarily map onto piecewise linear func- 
tions and thus the elements of the range are again approximated by the discrete 
operator 

K K 

F d u(x) ~J2Y1 F >v u jMx) - Fu(x) . (10) 

fc=l j=l 

In the following, for the sake of simplicity of notation, we identify the discrete 
operator F d with the matrix F of coefficients. Moreover, we assume that the 
discretization is fine enough that we can identify F and F d on PL(M). 

We use the following approximations for the fit-to-data term and the regu- 
larization functional: 

• Let the matrix A e M Arx3i be defined by the areas Ai of the triangles: 




if 3(i — 1) < j < 3i + 1, 
else. 
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Then, 



l\\Fu-y% HM) = \ JjFu(x)-y S (x)fdu g (x) 



(ii) 



t, jet, 
= i||A(VFu-V/)|| 2 

• In the applications presented below the rcgularization functional is either 
the total variation semi-norm or the quadratic Sobolev semi-norm of the 
gradient. We evaluate these functionals for u e PL(M) on the polyhedron: 

uePL(M)^ f \V M u\?dv(g), p = l,2. 

JM 

Let the diagonal matrix A £ M. LxL be defined by the areas A t of the 
triangles: 

Aii = 1Ai . 
From the above considerations we find that 



TZ(u) = - f \V M u\ p dv(g) 

P JM 

= -J22A i \(W M u) i \ p 

P T\ 

= -Tr(A(Z T Z)?) 
P 



(12) 



Because we have that 



— = VJ T G _1 JV. 
ou 



it follows that the derivative of the discrete functional T a y s (with M. replaced 
by M) at u is given by 

^j^(u) - ^F T V T A T A(VFu - Vy s ) 
ou 6 

+ a A(Z T Z)<"- 2 '/ 2 V T J T G- lT JV T (VJ T G iJVu) . 

The formal derivative of the dlZ(u) is a discrete approximation of the differential 
operator 

-div M (\V M u\ p ~ 2 V M u) ■ 
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In particular, for p = 2 we obtain a consistent approximation of the Laplace- 
Bcltrami operator. 

The optimality condition dT a y s (u, p) = for all p, can be solved with a 
Landweber fixed point iteration: 

u (*+i) = u (*) _ kV T q iV «(u<*>) k = 0,1,2,.... 

Here k denotes the step size and is chosen to satisfy a stability criterion [T5] . The 
algorithm is usually terminated if the difference of the update || u ( fc + 1 ) — u^ 5 )^ 
is below a given threshold for the first time. 



3.2 Applications 

Denoising of data on manifolds 

We consider denoising of image data on a closed finite dimensional Rieman- 
nian manifold. The usual assumption is that the data y s can be decomposed 
into a ideal image and additive white noise n a , with mean and variance a. 
That is 

V s = u* +n a with ||n CT || L =(M) < a ■ 

This corresponds to Equation (JlJ where the operator F is the identity. Thus 
denoising can be viewed as an inverse problem. The goal of denoising is to 
remove the noise component n a from y s but at the same time preserve the 
visual appearance of the clean image . 

In Figure [T] we compare quadratic Sobolev semi-norm regularization with 
total variation minimization. 



Image deblurring 

The general assumption is that the imaging data y s is obtained from the 
clean image by convolution with a smooth kernel function h and by additive 
white noise with mean zero and variance a. Thus, in the terminology of the 
paper, the operator equation reads as follows 

y d = h -k + n a =: Fu^ + n a . 

In our numerical experiments we assume that the kernel function h is a Gaufi- 
function with variance r. That is, 

1 fd 2 Jp,q)\ 
k ^ P ' ^ = 2^ 6XP I 2t 2 J f ° r eVery P,qe M ' 

where d g (p,q) denotes the geodesic distance on the polyhedron M. 



For implementing the Landweber algorithm (13) we use the discrete con- 
volution, which is (similar as in Section |3| written as a bold face matrix H 
with matrix entries d(vij,Vk,i)- The geodesic distance between two points p 
and q can be computed by solving the Eikonal equation with constant velocity 
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(a) u+ = Lenna on M (b) Testdata: lit + Af(0, a) on M 




(e) Result with total variation denoising (f) u s a — 



Figure 1: Denoising: The first row shows the ideal data u and the noisy signal 
y s . The second row depicts the results with total variation regularization and the 
residual image u s a — . The last rows results are obtained with by regularizing 
with the squared gradient. 
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Type SNR Original/Noise Original/Result 



Sobolev semi-norm 


Denoising 


22.12 


23.87 




Dcblurring 


16.89 


17.23 


TV semi-norm 


Denoising 


22.12 


26.35 




Deblurring 


16.89 


21.18 



Table 1: This table evaluates TV and quadratic variational regularization for 



denoising and deblurring in terms of the signal to noise ratio(c.f. 14 ) 



p(x) = 1 on M with an algorithm described in [T5]. That is, after fixing one 
point q e M, g?(x, q) solves the Eikonal equation: 

|Vd(i,g)| = p(x) • 

The Landweber algorithm for minimization of the discretized regularization 
functional reads as follows 

u (fc+i) _ u (*0 _ K (H T (Hu( fc > - y s ) + aV72.(u<*>)) . (13) 

Again, we compared TV and quadratic regularization. For TV regulariza- 
tion, k < 1+ a8/ e has to be chosen sufficiently small. Figure |2j shows results for 
deblurring with TV minimization and quadratic Tikhonov regularization. 

In Table[T]we summarized the results on the denoising and delurring problem 
for a fixed a. In order to compare the performance of different choices for 1Z we 
use the signal-to-noise ratio (SNR) measured in dB . The SNR is defined as 

SNR = 201og 10 f „ f t|L t 2 ,{ M) ). (14) 

The better performance of the Total Variation regularization stems from the 
fact that discontinuities along edges are preserved while the Sobolev semi-norm 
introduces severe blurring of the edges. 

Variational regularization for inversion of the spherical Funk-Radon 
transform 

In a recent work Louis et al [20] discuss a problem of density estimation, 
which requires the inversion of the Funk-Radon transform on the 2-sphere. In 
general, for arbitrary space dimension, the Funk-Radon transform maps a func- 
tion defined on the 2-sphere to its means over the great circles. That is, 

Fu(x) = — / u(y)dv(y) for every x £ § 2 . (15) 
w i Js 2 nx 1 - 
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(a) ut = Lenna on M (b) y s = Gaufiian convolution of + 

Af(0,<r) 




(c) Result with quadratic Sobolev semi- (d) it* — 

norm regularization 




(e) Result with total variation regulariza- (f) u s a — 

tion 



Figure 2: Deblurring: The first row shows the original signal and the noisy 
signal. The second row depicts the results with total variation regularization 
and the residual images, respectively. The last rows results are obtained with 
by regularizing with the squared gradient. 
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Helgason [T5] provids a closed form for the inverse of the Funk- Radon transform. 
In |20j an approximate inverse for the efficient numerical inversion of the Radon- 
Funk transform on the 2-sphere has been proposed. 

Here we investigate quadratic Tikhonov regularization with Sobolev semi- 
norm regularization term on the 2-sphere. The method consists in minimization 
of the functional 

T a , y >(u) = \\Fu - /||l2 (s2) +a|V Af w|^ ( g 2) . (16) 

The proposed numerical minimization algorithm requires real valued spherical 
harmonics : The functions Yj 1 "^,!/?), where I denotes the degree and m the 
order, form an orthonormal basis on S 2 : 

Y i m (x)Y e ™(x)dx= f f Yr(0,<p)YF'(O,<p)d(p,<p)=5u,6 mm ,. 
s 2 Je=oJ^=o 

In the following we define a single index j :— j(l,m) := (/ + 1)1 + m for I — 
0, 1, 2, . . . , L and m = 0, . . . , I and identify the coordinates x on the sphere with 
polar coordinates (9, tp). For the numerical minimization we use approximations 
of u 6 L 2 (M) with real spherical harmonics of maximal degree L. This is 

R 

u{9^)^Y, c i Y ^ 6 ^)- (17) 

j=i 

Therefore from discrete sample values u = (u^) = (it(xfc)) , k = 1, . . . , N x R 
the spherical harmonics expansions can be computed from the following matrix 
equation: 

u = Be , 



where 



B 



Y 1 (6 1 ,< Pl ) ... Y R {Q uVl ) 



Yi(9 N ,ip N ) ... Y R (0 N ,ip N ) 



is the matrix of spherical harmonics basis functions. The coefficients c are the 
coefficients of the best approximating solution in L 2 (M) and are given by 



c = (B B) B u . 

Figure [3] a) shows the best approximation of a function u with a spherical 
harmonics polynomial of degree 26. 

Using the Funk-Henke Theorem it has been shown in [7] that the Funk- 
Radon transform of a function u, given in a spherical harmonics basis, takes the 
simple form 

.R 

Fu(x)«27r£%)(0)c,-li(x). (18) 
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■ 



(a) Approximation of u' with spherical harmon- 
ics basis of maximal degree 26. 





(b) Approximation error. 




(c) y = Funk-Radon transform of (d) y s = y + additive Gaufiian noise 

Figure 3: Testdata and their Funk- Radon Transform 
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The Legendre polynomial of degree evaluated at is 



%)(0) 



if /(i) is odd, 

( _l ) i(i)/2 i-3-...( i (i)-i) ifZ(i)is even. 



2-4-6.. ./(j) 

Therefore, the discrete Funk-Radon transform can be written as 

F = BF , with Fjj = 2irP l{]) (0) 

a diagonal matrix. An example of the evaluation of the Funk-Radon transform 
can be seen in Figure [3] c). In Figure [3] d), the transformed signal is perturbed 
by some additive Gaufiian noise with variance er = 0.05. In order to reconstruct 
the signal from its Funk-Radon transform we minimize the Tikhonov functional 



from Equation ( 16 I 



As in the previous examples the subgradient of the L 2 (M)-norm of the gradi- 
ent the Laplace-Beltrami operator, which in spherical harmonics basis expansion 
is given by 

and in matrix notation 

V j = + I)- 

Reconstruction of the inverse Funk-Radon transform with Tikhonov regulariza- 
tion requires solving the linear system 

(F T B T BF + aL)c = F T B T y 5 . (19) 

This equation can be solved again with a Landweber iteration. This equation 
has been solved again via a Landweber iteration: 

g(n+l) = £ (n) _ K ( F T B T (BF£ („) _ y 5) + aL g(n^ (20) 

Since the Funk transform annihilates odd functions (see |18|), we take an 
even function to test our inversion algorithm. As in |20) we use the function 

u(x) = cos(37r(z — y)) + cos(37rx); 

and evaluate the function at 900 point on the sphere as they are provides in 
[25] . In the numerical experiments we used spherical harmonics of degree 26. 
The reconstructions are depicted in Figure [4] Note that if a = we cannot 
solve for C because the matrix F does not have full rank. Only by regularizing 
the inversion with the Laplace-Beltrami operator allows the reconstruction of 
u(x). In Figure [1] we can observe the smoothing effect of the Laplace-Beltrami 
operator on the solution. In the left column of Figure|4]a low value a = 0.3 and 
in the right column the result with a high value a = 1.2. We observe a much 
smoother reconstruction for higher values of a. 
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(a) Reconstruction with low regularization 



(b) Reconstruction with strong regularization 
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4 Conclusion 



In this paper, we have studied the problem of variational regularization of 
inverse and ill-posed problems for functions on closed Riemannian manifolds. 
The analysis (stability, convergence, and rates) follows from standard results 
on convex regularization and are reviewed. The main contribution of this pa- 
per concerns the numerical analysis of such regularization methods and the 
numerical implementation. Moreover, three inverse problems appearing in non- 
destructive evaluation and Computer Vision are discussed. 
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5 Background and definitions 

All along this paper we use the notation: 

• Ai denotes a closed n-dimensional manifold in 

• M = (V, T) denotes a polygonal approximation of Ai. In this paper, 
typically, it is a polyhedron. 

• x denotes coordinates on the manifold, x denotes coordinates in the Eu- 
clidean space. 

• d denotes the subdifferential, V is the gradient in the Euclidean setting, 
Vm an d V^vi are the covariant derivatives on the Riemannian manifolds 
M, Ai, respectively. 

• If not specified otherwise |.| denotes an Euclidean distance. 

In the following we review elementary facts from Riemannian geometry and 
nonlinear analysis on manifolds. 

1. The metric tensor gij{x) expressed in a coordinate chart (f2, ip) is 

°-*M-<^^>«- ■ 

The inverse metric tensor is g % i = G _1 . 

2. Given a smooth, closed Riemannian ra-manifold (Ai,g), there is an as- 
sociated positive Radon measure on Ai, the Riemannian measure, which 
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is defined as follows: Given an integrable function u : M. h-> K, an atlas 
(f2j, tfi)i<£i of .A/f and a partition of unity (fij, t/?j, aj)j e j, 



/ udv(g) = ^2 (ajuo <p jy /\f\)dx 



(21) 



where = |det(<?.y)| is the volume of the metric tensor and dx is the 
Lebesgue volume element on R" and therefore dv(g) = y / Jg\dx 



3. The gradient can be expressed with the chart (Q, ip) 

d(u o ip^ 1 ) 



(V M u)i = g v] 



dxj 



(22) 



4. The adjoint operator of the gradient is the divergence, which satisfies for 
a given vector field X on A4: 



/ udiv_wXdv(g) = — V m u ■ Xdv(g) 
J M Jm 



(23) 



In a chart (fi, ip) the divergence is obtained by using Equation 22 in Equa- 
tion [231 



divjviX 



1 



i=l 



(24) 



5. Given (A4, g) and 7 : [a, 6] i-> M a curve on j\4, then £(7) is the length of 
the curve on M. with respect to g. For p, g on X with 7(a) = p 7(a) = p 
and 7(6) = q, the distance associated with g between two points p and q 
is 

d g (p,q) = inf £(7). 

76 c pq 

C pq is the space of continuous curves connecting p and q. The distance 
d g defines a metric on the manifold. In this paper, we assume that the 
metric space (M,d g ) is always complete. 

6. In a closed Riemannian manifold (without boundary) the Hopf-Rinow 
theorem implies that for every pair of points on the manifold there exists 
a unique geodesic [3]. 

Given a smooth and closed n-dimensional Riemannian manifold (AA,g), we 
define (see [17] ) 



C{{M) :-- 



eC°°(M):\\u\\ k , p := \\u\\ p + £ 



.VI 



d%) 



< 00 
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where |.| p is the L p (A^)-norm with respect to the Riemannian measure dv(g). 
The space W p ' k (M) (see e.g. [T7]) is defined as the completion of the space 
C%(M) with respect to the norm \\.\k,p- In particular for k = 1 and p > 1, we 
have 

H|vkp>i(A4) = |«| p + |V^w|p ■ 

Now, we recall the definition of the space of functions of bounded variation on 
manifolds. 

Definition 5.1. We define BV(A4) as the space of functions with bounded 
variation and is the set of functions u G L 1 (A / J) such that \Dm u \ (M) < +oo. 
The space is endowed with the norm \u\gv{M) = IMIl 1 ^) ~H-^jV( u I > where 
\Dj^u\ (M) denotes the variation of u, which is defined by Q) 

The space BV(M) is a Banach space endowed with the norm ||.||_bv(a-i)- 
It can be understood as the natural (weak) closure of W 1:1 (A4). Due to the 
theorem of Meyer and Serrin |21j it is possible to approximate Sobolev functions 
defined on subset of the Euclidean space by smooth functions. For the sake of 
completeness we provide a proof to show the essential difference in the manifold 
setting. 

Theorem 5.2 (Approximation of BV- Functions). Let M. be a smooth, closed 
Riemannian manifold and u £ BV{M). Then there exists a sequence (u n ) of 
functions in C^°(M) such that 

• u n — > u in l}(M), 

• Im |Vjwu„(x)| dx \D M u\ (M). 

Proof. The proof is closely related to f^j, where weighted BV spaces have been 
considered and thus omitted. 

The second important property of BV functions used in this paper is covered 
by the following embedding theorem: 

Theorem 5.3 (Compactness Theorem). Let M be a closed manifold, and let 
(ti„)„ be a sequence of functions inBV(M) such that sup„ \Dj^u n \ (M.) < +oo. 
Then there exists a subsequence of u G BV(A4) converging strongly in L 1 (A4). 

Follows from combining the analogous result for functions W 1,1 (A4), which 



is stated in HebevfT?]. and Theorem 5.2 



Theorem 5.4 (Embedding Theorem). For every function u € BV(Ai) 

M h ^r (M) <C(n)\D M u\(M). (25) 
The proof is analogous to the Euclidean setting and thus omitted. 
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